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Kinetics of crystal-growth is investigated along the solid-liquid coexistence line for the (100), 

(110) and (111) orientations of the Lennard-Jones and Weeks-Chandler-Andersen fee crystal-liquid 
interface, using non-equilibrium molecular dynamics simulations. A slowing down of the growth 
kinetics along the coexistence line is observed, which is mostly a temperature effect, with other 
quantities such as the melting pressure and liquid self-diffusion coefficient having a negligible impact. 

The growth kinetics of the two potentials become similar at large values of the melting temperature 
and pressure, when both resemble a purely repulsive soft-sphere potential. Classical models of 
crystallization from the melt are in reasonable qualitative agreement with our simulation data. 

Finally, several one-phase empirical melting/freezing rules are studied with respect to their validity 
along the coexistence line. 


I. INTRODUCTION 

A central quantity for the understanding of crystalliza¬ 
tion processes from the undercooled liquid is the kinetic 
growth coefficient [I|. The kinetic growth coefficient, 
is defined as the constant of proportionality between the 
velocity Vi with which the crystal-liquid interface moves 
and the interfacial undercooling, AT = Tm — T, 

Vi = pAT, ( 1 ) 

with Tm the melting temperature. Note that Eq. o is 
only expected to hold if the undercooling AT is suffi¬ 
ciently small. Magnitude and anisotropy of the kinetic 
growth coefficient play a dominant role in determining 
the morphology of the growing crystal ii and are also 
essential parameters required for the continuum mod¬ 
elling of solidification processes Q. 

Experimental measurements of the kinetic growth co¬ 
efficient have been scarce, with the exception of a few 
studies on metallic systems Si and white phosphorous 
(P4) 0. However, atomistic simulation techniques 0,0 
such as Molecular Dynamics (MD) provide detailed in¬ 
formation about the microscopic structure and dynamics 
of the interface region. Thus, atomistic simulations can 
be used to test various analytical approaches to describe 
crystal growth such as the Wilson-Frenkel I, 0 and 
Broughton-Gilmer-Jackson PS model. 

Different simulation techniques have been developed 
for the investigation of crystal growth kinetics |lll - l^ . 
In the capillary fluctuation method (isl - f^ . the kinetic 
growth coefficient /i is obtained from equilibrium MD 
simulations by analyzing the height fluctuations of the 
crystal-liquid interface at coexistence. Using this ap¬ 
proach, /i has been computed for hard spheres [10 , met¬ 
als El, a Lennard-Jones system and the TIP4P/2005 
water model [20. A variation of the capillary fluctuation 
method has been proposed by Tepper and Briels [20-[20 . 
In their approach, /r is extracted from the equilibrium 
fluctuations of the number of crystalline atoms in an in¬ 
homogeneous solid-liquid system. 


Another widely used approach to obtain kinetic growth 
coefficients is the free solidification method (FSM) P, 
[20 which is based on non-equilibrium MD. Here, one sim¬ 
ulates inhomogeneous systems where the crystal is sepa¬ 
rated from the liquid phase via two planar interfaces (two 
interfaces appear due to periodic boundary conditions). 
By monitoring the rate of change of the system’s volume 
with respect to time, the kinetic growth coefficient, /i, can 
be determined. This approach has been applied to var¬ 
ious one- and two-component metals El [2a - [30 as well 
as Lennard-Jones systems [2l,[30|. The estimated values 
of fjL for metallic systems, obtained from FSM, have been 
shown to be in good agreement with those obtained for 
hard spheres from the capillary fluctuation method [10 . 

For systems with a crystal face-centered cubic (fee) 
phase, it has been suggested from the latter simula¬ 
tion studies that the diffusion-limited classical Wilson- 
Frenkel model of crystal growth has to be modified into 
a collision-limited growth model to explain the high crys¬ 
tal growth rates corresponding to the (100) orientation 
of the crystal-liquid interface. Moreover, the collision- 
limited growth model seems to be a good predictor of 
for the (110) interface of metallic systems [ 10 ], too. Only 
the (111) interface tends to follow the Wilson-Frenkel ki¬ 
netics (l0 . 

Most of the above studies have been done under “am¬ 
bient conditions”, i.e. at the melting temperature cor¬ 
responding to zero pressure conditions. Little is known 
about the dependence of the kinetic growth coefficient on 
pressure and temperature along the coexistence line. In 
this work, we investigate the growth kinetics of the (100), 
(110), and (111) orientations of a fee crystal-liquid inter¬ 
face for two different models, employing the FSM: (i) 
a force-shifted Lennard-Jones (fsLJ) potential and (ii) a 
purely repulsive Weeks-Chandler-Andersen (WCA) po¬ 
tential. For both models, systems at various pressures 
and temperatures along the coexistence line are studied. 
From the FSM simulations, the coexisting temperatures 
and pressures as well as the kinetic growth coefficients 
are obtained. 
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Our results indicate that the growth kinetics depends 
only weakly on pressure while it is significantly affected 
by a change of the melting temperature. For the fsLJ po¬ 
tential, there is an initial regime where the coexistence 
pressure changes by two orders of magnitude while the 
melting temperature remains essentially unchanged. In 
this regime, the kinetic growth coefficients are almost 
constant as a function of pressure. For both the fsLJ 
and the WCA models, however, an increasing melting 
temperature Tm (with increasing the “melting” pressure 
Pm) is associated with a slowing-down of the growth ki¬ 
netics. At high values of Pm and Pm, the kinetic growth 
coefficients in reduced units tend to reach asymptotic val¬ 
ues which correspond to the ones found for hard spheres, 
though not identical. In our analysis, we discuss to 
what extent the collision-limited crystal growth model 
of Broughton, Gilmer, and Jackson (BGJ) is valid and 
relate the growth kinetics to the self-diffusion coefficient 
of the liquid, the entropy of fusion and the liquid and 
crystal coexistence densities. In this context, we also 
discuss empirical rules of melting and freezing [H, [s^ , 
namely the Lindemann [s^, the Raveche, Mountain and 
Streett and the Hansen-Verlet criterion. 

In the next Section, we describe the two interaction 
potentials considered in this work. In Section III, we 
outline the FSM, followed by Section IV on the simula¬ 
tion details. The results are presented in Section V, and 
finally we come to the conclusions in Section VI. 


II. INTERACTION POTENTIALS 


Simulations are carried our for two different models 
that are derived from the Lennard-Jones (LJ) potential. 
This potential describes the interaction between two par¬ 
ticles separated by a distance r by the function 


^(r) = 4e 



with e and a being two parameters setting respectively 
the microscopic energy and length scales for two neigh¬ 
boring particles. 

The first model considered in this work is a force- 
shifted Lennard-Jones (fsLJ) model, defined by 


UtsLj{r) 


(j){r) - (j){rc) - 4''{rc)[r - rc] r < 

0 otherwise 


( 3 ) 

with (j)' = d(j)/dr. The cut-off of the potential is set to 


Tc = 2.5 a. 

The second model is a Weeks-Chandler-Andersen 
(WCA) potential, i.e. a LJ potential which is cut off at 
its minimum at = 2^/® cr and shifted to zero. 


UwCA{r) 


4>{r) + e r < 2^^® a 

0 otherwise. 


( 4 ) 


Thus, the WCA model is a purely repulsive potential. In 
the fsLJ model, only at very high coexistence pressures 


the repulsive part is expected to dominate kinetic proper¬ 
ties and phase behavior and so we can study how attrac¬ 
tions between the particles affect crystal growth along 
the coexistence line going from low to high coexistence 
pressures. 

In the following, energies and lengths are expressed in 
units of e and cr, respectively, and the masses of the par¬ 
ticles are set to to = I. Thus, thermal energy, k^T (with 
the Boltzmann constant fee = 1 and T the temperature), 
and pressure, P, are expressed in units of e and e/tr®, 
respectively. Time is measured in units of r = 
while the kinetic growth coefficient is reported in units 
of k^/yjme. 


III. FREE SOLIDIFICATION METHOD (FSM) 

A crystal in contact with its melt at a temperature T 
below (above) the melting temperature Tm (cf. Fig. [T]) 
will grow (shrink) until the entire system crystallizes 
(melts). The FSM [II, [H, [2^ [sl, Q is an approach 
to compute the crystal growth (or melting) rate as well 
as the coexistence temperature Tm at a given pressure P. 
The starting point are standard isothermal-isobaric MD 
simulations at various temperatures T and at a particular 
value of the pressure, P, keeping the number of particles 
N constant. From these NPT simulations, the tempera¬ 
ture dependence of the density of the crystal (in our case 
a fee crystal) as well as the melt are determined. Grad¬ 
ually increasing the temperature of the crystal leads to 
melting at a temperature Ti while the subsequent gradual 
reduction of the temperature leads to re-crystallization at 
a different temperature T 2 . Thus, hysteresis is observed, 
i.e. the heating and cooling curves do not follow the same 
path. Such “heating-cooling” plots indicate the region in 
which the melting temperature, Tm, is located. 

Now, the FSM scheme consists of the following 
steps m, First, N atoms are arranged on a fee 

lattice in an elongated simulation box of size L^xLyX 
with the desired orientation of the crystal pointing in z 
direction and the length of the system along the z direc¬ 
tion being approximately five times that in x and y direc¬ 
tions (cf. Fig. [ 1 ]). At each temperature and pressure, the 
density of the fee crystal is obtained from the aforemen¬ 
tioned heating-cooling plots (Fig. [2|). Then, at a given 
temperature and pressure the system is equilibrated in a 
NPxPyPzT ensemble [s^ (with T^, Py, and Pz respec¬ 
tively the pressures along the a;, y and z directions) in the 
range in which hysteresis is observed. The reason for car¬ 
rying out simulations in the constant pressure ensemble 
is to ensure that the crystal is free of any residual stress 
along the three Cartesian axes. Moreover, for maintain¬ 
ing constant pressure along the different Cartesian axes, 
simulations are carried out in the NP^PyPzT ensemble 
rather than in the NPT ensemble since the simulation 
box is a cuboid with unequal lengths along the different 
Cartesian axes. After equilibration is reached, the aver¬ 
age length of the simulation box in the x and y directions 
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are determined. 

After relaxing the crystal sample in the first step, two- 
thirds of the atoms in the middle of the box are fixed 
and the rest of the system is heated up to a high tem¬ 
perature, T » Tmi to eventually melt it. In this step, 
the lengths of the simulation box in the x and y direc¬ 
tions are fixed to the average lengths obtained from the 
previous NP^PyP^T ensemble run, and the simulation is 
carried out in the NPzAT ensemble by varying only the 
length Lz (i.e., maintaining constant pressure Pz). Here, 
A corresponds to the area of the system {A = x Ly). 
In the third step, the temperature of the whole system is 
set back to the initial temperature in which the crystal 
was prepared, with the atoms in the middle of the region 
still fixed. This simulation in the NPzAT ensemble runs 
for a short period, just long enough to cool the melted re¬ 
gion to the desired temperature. Finally, all the particles 
are allowed to move and the simulation in the NPzAT 
ensemble is continued. 

In the steady-state, the length of the system, varies 
linearly with time and one obtains the change of Lz per 
unit time, Lz, from the slope of Lz(t). From Lz and a 
mass balance equation, it is straightforward to obtain the 
interface velocity via [35l| 



t=0 







t=2500 


P\Lz 

2(Pc - Pi) 


( 5 ) 


with Pc and p\ the bulk densities of the solid and liquid 
phase, respectively. Close to the melting temperature 
the interface velocity varies linearly with temperature. 
From a linear fit to the temperature dependence of the 
interface velocity, the kinetic growth coefficient p as well 
as the melting temperature Tm are obtained, cf. Eq. ©• 


IV. SIMULATION DETAILS 

To integrate the equations of motion, the velocity Ver- 
let algorithm Q is used. Pressure is kept constant by 
coupling the system to an Andersen barostat [3, |4l| . 

The coupling parameter of the barostat (the mass of the 
piston) is set to M = 0.001 for the NPT runs to deter¬ 
mine the “heating-cooling” curves and to M = 10 for all 
the FSM simulations in the NP^PyP^T and the NP^AT 
ensemble. To keep the temperature constant, the system 
is coupled to a stochastic heat bath by assigning every 
200 time steps random velocities to the atoms sampled 
from a Maxwell-Boltzmann distribution. The reduced 
time step is taken to be 5t = 0.005r for all the simula¬ 
tions. 

To compute the density of the crystalline and liq¬ 
uid phases at coexistence and to determine the density- 
temperature hysteresis curves, NPT simulations with 
N = 2048 particles are carried out at various values of 
the pressure in the interval 0.005 < P < 32 for the fsLJ 
potential and in the interval 1.0 < P < 1020 for the 
WCA potential. Initially, the particles are placed on a 


FIG. 1. (Color online) Time evolution of a crystal-liquid inter¬ 
face during growth (T < Tm), corresponding to the (100) ori¬ 
entation of the fsLJ system at the pressure P = 1.0 and tem¬ 
perature T — 0.65, represented by snapshots of the system at 
various times t (the melting temperature is Tm = 0.653). Ini¬ 
tially (topmost configuration), there is a crystal sandwiched 
by equal amounts of liquid on both sides. With time the liq¬ 
uid portion shrinks gradually and eventually the entire system 
crystallizes. Time t is measured in units of r (see text). 



FIG. 2. (Color online) Heating-cooling curves in the density- 
temperature plane corresponding to (a) the fsLJ and (c) the 
WCA model at various applied pressures. Also plotted is the 
temperature dependence of for the fsLJ and the WCA 

model in panels (b) and (d), respectively. The dotted and 
dashed lines in panels (b) and (d) represent the coexistence 
values of the crystal and liquid, respectively. 
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fee lattiee in a eubie eell of dimensions Lx Lx L. We equi¬ 
librate the system for 25000 time steps and then perform 
produetion runs for another 25000 steps. From the data 
eolleeted during the produetion runs, the volume of the 
system is determined to obtain the equilibrium density. 
The temperature of the system is raised by a small step 
and the above proeedure is repeated to obtain the density 
for the next higher temperature. This process is contin¬ 
ued until the crystal melts to form the liquid phase. The 
same procedure is followed to obtain the bulk liquid equi¬ 
librium density by gradually lowering the temperature of 
the system in small steps and calculating the equilibrium 
density at each temperature. 

For the FSM simulations of the fsLJ model, parti¬ 
cles are placed in an elongated box of size L x L x 5L. 
Systems containing N = 34560, 35088 and 35700 parti¬ 
cles are considered respectively for the (100), (110), and 
(111) orientation of the fee lattice in z direction. With 
the same relative lengths of the simulation box along 
the X, y and z directions, simulations with N = 14580, 
15444 and 15795 particles are carried out for the WCA 
model along the (100), (110), and (111) orientation, re¬ 
spectively. At several temperatures in the hysteresis re¬ 
gion (see Sec. IIIIll crystals are equilibrated for 75000 
time-steps in the NPy^PyP^T ensemble. From the last 
10000 time-steps, the average lengths of the simulation 
cell along the x and y directions is determined. Ly and 
Ly are fixed to this value to carry out simulations for the 
next step in the NP^AT ensemble. 

In the second step, the crystalline particles in the mid¬ 
dle one-third of the box are fixed, while particles in the 
remaining region are equilibrated at a high temperature 
for 150000 time steps to melt the system. Then, keeping 
the crystalline region fixed, the liquid is cooled to the 
desired temperature in a short run of 10000 time steps. 
Finally, all the particles are allowed to move, performing 
a run over 500000 time steps from which the interface 
velocity Vi is determined. Statistical errors were deter¬ 
mined from 10 — 20 independent realizations. 


V. RESULTS 

A. Heating-cooling plots 

FiguresOi and c show heating-cooling plots of the den¬ 
sity at different pressures for the fsLJ and the WCA 
model, respectively. With increasing pressure, hystere¬ 
sis is observed in a larger temperature interval. As we 
shall see below (see Sect. IVEl) , this is associated with 
a gradual slowing down of crystal growth kinetics. By 
scaling the density by a factor 1/T^/^, all the curves cor¬ 
responding to the various pressures tend to fall in a sim¬ 
ilar range along the j/-axis (see Fig. and d). Also, the 
coexistence values of the quantity for the crystal 

and liquid asymptotically reach a similar constant value 
for both the fsLJ and WCA model, indicating that the 
phase behavior of both systems at high temperature and 


pressure is similar. We discuss this point in more detail 
in Sect. lYSl 


B. Interface velocity 

The behavior of the crystal-melt interface depends on 
the temperature at which it is simulated. The crystal will 
grow below the melting temperature while above it, melt¬ 
ing will occur. Since the crystal density is higher than 
the melt, the length of the system along which the inter¬ 
face is oriented will increase during melting and shrink 
during crystallization. In Fig. [3^, we plot the length of 
the system along the z direction, versus time for dif¬ 
ferent temperatures corresponding to the fsLJ potential 
at P = I.O. The data is averaged over 10 independent 
realizations. Figure shows that after a transient pe¬ 
riod at the beginning, varies linearly with time when 
the steady state is reached and ultimately reaches a con¬ 
stant at long times, when the whole system has either 
crystallized or melted. Just prior to this, we find a non¬ 
linear regime where the crystallization and melting are 
much faster than in the linear steady-state regime, be¬ 
cause one of the phases has shrunk to such a small size 
that the two interfaces interact with each other (see the 
two bottom-most curves in Fig. |3^). 

Figure[3]3 shows v\ as a function of temperature for the 
three different crystal orientations i.e. (100), (110) and 
(111) at P = 1.0. The melting temperature is dictated 
by thermodynamics and is expected to be identical for all 
crystal orientations. At small undercoolings, the system 
is in the linear response regime. Hence, around Tm, the 
simulation data for Ui can be fitted by a linear law m = 
/^(T’m — T), with Ui = 0 at Tm- From Fig. we find 
that for the (100) crystal orientation, Vi vanishes around 
T = 0.653. For the (110) and (111) crystal orientations, 
Ui approaches zero respectively at T = 0.653 and T = 
0.655. 


C. Thermodynamic properties at coexistence 

The FSM simulations were carried out at different pres¬ 
sures to obtain the respective melting temperatures. Fig¬ 
ure m shows the phase diagram of the fsLJ system and 
the WCA potential along the P — T plane. Our simu¬ 
lation data corresponding to the coexistence conditions 
is in very good agreement with that obtained previously 
using Gibbs-Duhem integration from a known melting 
temperature and pressure at a single coexistence point 
for the fsLJ and the WCA model [i^. At low pres¬ 
sures, there is a significant difference between the coex¬ 
istence curves of the two potentials on account of the 
different roles played by the — 1/r® term. Due to this 
attractive part of the fsLJ potential, the particles sit at 
the potential well at low pressures (Pm ^ 0.5) and as 
a result the melting temperature stays almost the same 
even when the potential changes by two orders of mag- 
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FIG. 3. (Color online) (a) Length of the system along the 2 : 
direction, Lz vs. time for different temperatures at P = 1.0. 
(b) Interfacial velocity as a function of the temperature at 
P = 1.0 for the (100) (circles), (110) (squares), and (111) 
(diamonds) crystal orientations. Straight lines are linear fits 
to the data at low undercoolings (see text). 



FIG. 4. (Color online) Phase diagram in the pressure- 
temperature plane for the fsLJ and the WCA model, as indi¬ 
cated.The dotted and dash-dotted lines represent coexistence 
values obtained by Ahmed and Sadus for the WCA poten¬ 
tial and by Errington et al. for the fsLJ potential fi^ . 
The inset shows the quantity as a function of the 

melting temperature Tm. The green solid line corresponds to 
the melting pressure for the WCA potential at Tm = 1.0. 
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FIG. 5. (Color online) Variation of f.d.c. and r.d.d. (inset) as 
a function of the coexistence pressure for the fsLJ and WCA 
potential. The green dashed horizontal lines represent the 
corresponding values for the inverse twelfth-power potential 
(see text). 


nitude (from P = 0.005 to P = 0.5). For Pm > 0.5, the 
melting temperature increases rapidly with the pressure 
as the repulsive part of the potential becomes dominant 
while the attractive term plays less and less of a role in 
determining the phase behavior. One can clearly identify 
these two regimes in the coexistence line corresponding 
to the fsLJ potential as shown in Fig. 0] 

Figured] also shows that the P — T coexistence curve 
of the WCA model is almost a straight continuous line 
in the whole considered range of melting temperatures. 
Moreover, at high values of Tm the coexistence line of 
the fsLJ model seems to become identical to that of the 
WCA model. This behavior is expected because at high 
temperatures the phase behavior of both models is dom¬ 
inated by the repulsive 1/r^^ interactions, in agreement 
with the findings in Refs, [dsl - ld^ . 

It is interesting, therefore, to compare the coexis¬ 
tence behavior of the fsLJ and the WCA potentials with 
the inverse 12th-power soft-sphere potential, U{r) = 
e{alr)'^ {n = 12). Inverse-power law potentials are 
fully determined by one parameter (here ecr") and co¬ 
existence is fully specified by a single quantity, r„ = 
As a consequence, the reduced melt¬ 
ing pressure shows a power-law scaling with respect to 
the reduced melting temperature. Pm = PiTm^'*’^'^”. For 
n = 12, this relation reduces to Pm = PiTm^'^^, where Pi 
is the pressure corresponding to Tm = 1.0. In the inset 
of Fig. dl we plot the quantity Pm/Tm^'^^ as a function 
of the melting temperature for both the WCA and the 
fsLJ model. In case of the WCA potential (Pi = 13.41), 
this relation is satisfied for Tm > 1.0. For the fsLJ 
potential, at larger melting temperatures, this quantity 
tends to the same value as in case of the WCA potential 
i.e. Pm/T]Jj-^® = 13.41. However, the melting tempera¬ 
ture corresponding to this value of the pressure is around 
Tm = 1.48, indicating that for Tm > 1.48 (or Pm > 14.0) 
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FIG. 6. (Color online) Entropy of fusion as a function of 
the melting pressure corresponding to the WCA and the fsLJ 
system. The green dashed horizontal line represents AS' for 
the inverse twelfth-power potential (see text). 


the phase behavior of the fsLJ potential approaches that 
of the purely repulsive inverse 12th-power soft-sphere po¬ 
tential. 

In Fig. [5)3 and d, we have plotted the scaled density 
as a function of temperature that tends to ap¬ 
proach a constant value in the high-temperature limit 
that corresponds to that of a 1/soft-sphere potential. 
For the latter potential, the estimate of Hoover et al. [i^ 
for the crystal (liquid) coexistence value of is 

0.844 (0.813), while for the fsLJ and WCA interaction 
potentials the coexistence values of the crystal (liquid) 
at the highest pressures considered are 0.941 (0.889) and 
0.873 (0.838) at Pm = 32 and 1020, respectively. Thus, 
also with respect to the scaled density both the 

fsLJ and the WCA model approach the value obtained 
for the soft-sphere potential. 

Along the coexistence line, we now discuss further 
thermodynamic properties, namely the fractional density 
change at freezing, f.d.c. = {pc — p\)/p\, also known as the 
miscibility gap, the relative density difference at freezing, 
r.d.d. = 2(pc — Pi)/(Pc + Pi), and the entropy of fusion 
AS = I/Tm- Here, I is the latent heat as obtained from 
the difference of the enthalpies of liquid and crystal at 
coexistence. 

In Fig. (5) we show the ratios f.d.c. and r.d.d. (inset) 
as a function of the melting pressure for the two interac¬ 
tion potentials. Both the miscibility gap and the relative 
density difference at freezing decrease along the coexis¬ 
tence line. Data corresponding to the WCA potential 
are in good qualitative agreement with those obtained in 
an earlier work [4^ but the magnitudes of the f.d.c. and 
r.d.d. reported by us are slightly higher. At larger pres¬ 
sures the fsLJ data tends to attain similar values as those 
of the WCA potential though somewhat lower. The 
f.d.c. and r.d.d. ratios corresponding to the poten¬ 
tial are 0.038 and 0.037 [i^, respectively (indicated by 


FIG. 7. (Color online) The self-diffusion coefficients D (main 
figure) and D^* (in the inset) as a function of the coexistence 
pressure Pm corresponding to the fsLJ potential. 

the horizontal lines in Fig. [5]). Figure El clearly shows that 
at large pressures both the WCA and fsLJ data converge 
to the same values as those of the r~^‘^ potential. 

The entropy of fusion (see Fig. E) decreases with 
increasing temperature and pressure for both interac¬ 
tion potentials in qualitative agreement with a previous 
work [i^. This indicates that at lower melting temper¬ 
atures there is a greater positional ordering of the solid 
phase as compared to that at higher temperatures. At 
the lower pressures, the entropy of fusion for the fsLJ po¬ 
tential approaches a constant value as the melting tem¬ 
perature changes little. At larger pressures, the values of 
AS converge to that of the potential, AS = 0.9 [4^ . 

We have seen that with respect to the thermodynamic 
properties, the fsLJ and the WCA model become simi¬ 
lar in the high-temperature limit where, in both cases, 
repulsive interactions dominate the phase behavior. 
Below, we show (Sect. FV^ that the crystal growth kinet¬ 
ics of the fsLJ and the WCA model also becomes similar 
at high temperatures/pressures. 

D. Self-diffusion coefficient of the liquid 

One of the classical models of crystal growth kinet¬ 
ics is the Wilson-Frenkel model [^, [l^ which describes 
crystal growth by an activated process, limited by the 
self-diffusion of the atoms in the liquid phase. Thus, this 
model predicts that the diffusion dynamics of the liquid 
strongly affects the growth kinetics. Therefore, we now 
analyze the self-diffusion coefficient of the liquid phase for 
temperatures and pressures along the coexistence line. 

The self-diffusion coefficient, D is computed from 
the mean-squared displacement of a tagged particle Q, 
6r'^{t) = ((rtag(t)-rtag(O))^) (with rtag(t) the position of 
the tagged particle at time t), via the Einstein relation, 
D = limt^oo Sr^{t)/6t. 

In Figs. |7| and E) we display D along the P — T coexis- 
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FIG. 8. (Color online) The self-diffusion coefficients D (main 
figure) and D* (in the inset) as a function of the coexistence 
pressure Pm corresponding to the WCA potential. 



fence line for the fsLJ and the WCA model, respectively. 
As is evident from Fig. [71 D for the fsLJ model remains 
unchanged in the pressure range 0.005 < Pm < 0.5, 
where the corresponding melting temperature changes 
very little. However, when the melting temperature sig¬ 
nificantly increases by about a factor of three the dif¬ 
fusion coefficient also increases by about 50%. For the 
WCA potential (Fig.|8|), the diffusion constant increases 
with increasing pressure due to the change in Pm all along 
the coexistence line. Overall, the increase in diffusion 
coefficient with increasing melting pressures (provided 
Tm increases), might indicate that the growth kinetics 
become faster along the coexistence line. However, as 
reported below in Sect. IV El the kinetic growth coef¬ 
ficient decreases at high pressures, indicating that the 
self-diffusion coefficient of the liquid does not play a sig¬ 
nificant role in determining the growth kinetics. 

For comparing the fsLJ and WCA potentials, the 
scaled diffusion constant, D* = D/Dsc (with Dsc = 
being a natural MD scale), is shown in 
the insets of Figs. |7|and|51 D* for the fsLJ potential de¬ 
creases with increasing pressure and approaches a value 
between 0.029 and 0.030 at large pressures. In case of the 
WCA potential, D* slightly increases at low pressures 
and then, in the pressure range 5.5 < P < 1020, remains 
almost unchanged, also at a value between 0.029 and 
0.030. Thus, also with respect to the scaled self-diffusion 
coefficient D* similar values are obtained for the fsLJ and 
the WCA model at high temperatures/pressures. 


E. Crystal growth kinetics 

The kinetic growth coefficient p is extracted from the 
slope of the linear fit to the interfacial velocity Vi at small 
undercoolings. Our results for both the fsLJ and WCA 
model confirm observations of prior studies [l^ [l9l| re¬ 


garding the magnitude of p for the different orientations: 
/^loo > Miio > Pill- The data for the three different crys¬ 
tal orientations along the coexistence line are reported in 
Figs. and [TUh for the fsLJ and WCA model, respec¬ 
tively. For the fsLJ model, p remains essentially constant 
in the pressure range 0.005 < P < 1.0 where the melting 
temperature changes only weakly. For larger pressures, 
when the melting temperature increases rather sharply, 
the kinetic growth coefficient decreases, indicating this 
decrease to be of a thermal origin. Similarly, in case of 
the WCA potential, p decreases with increasing pressure 
reflecting the variation of the melting temperature with 
respect to the coexistence pressure (Fig.lTOh). At a pres¬ 
sure of Pm ~ 30, the values of p for both models are 
already in good quantitative agreement. 

To compare the values of p for the fsLJ and the WCA 
model with those of the hard sphere system, we have also 
computed the reduced coefficient, p* = p^/ksTM/m or, 
with m = = 1, p* = p^/Tm- It is to be noted that the 

phase behavior of hard spheres is purely entropy-driven 
and interfacial properties are solely determined by pack¬ 
ing effects. In Figs. |9]3 and [TOb . p* as a function of the 
melting temperature is shown for the fsLJ and the WCA 
model, respectively. In agreement with a recent simula¬ 
tion study [i^, for both models p* tends to approach a 
constant value in the high temperature limit. The hard- 
sphere values of p* , as obtained by Amini and Laird [a, 
are pioo = 1.44(7), pno = 1.10(5) and pm = 0.64(4) for 



FIG. 9. (Color online) The kinetic growth coefficient p cor¬ 
responding to the fsLJ potential for the three crystal orien¬ 
tations (100) (circles), (110) (squares), and (111) (diamonds) 

(a) as a function of the coexistence pressure Pm and (b) as 
a function of the coexistence temperature. In (a), the solid 
lines represent the fits obtained from the classical models. In 

(b) , p* — py/Tm is shown. 


















the (100), (110), and (111) orientation, respectively. The 
asymptotic values of the fsL J and the WCA model in the 
high pressure/high temperature limit are close to these 
values, though not identical. 

The ratios /iioo/Miio and /rioo//.tiii reported in 
Ref. [l^ for hard spheres are 1.31 ± 0.09 and 2.25 ± 0.18, 
respectively. For a LJ system, a previous simulation 
study [i^ reports the values 1.53 for /iioo/Mno and 1.996 
for //ioo/mih at P = 0.00254. In Fig. [TTl we plot 
Mioo//^iio and /j.ioo/miii as a function of Pm for the fsLJ 
model and in the inset for the WCA model. For the 
fsLJ model, we find that both ratios are lower than the 
hard-sphere and LJ values in the pressure range where 
the melting temperature hardly changes. At higher val¬ 
ues of the melting pressure, both ratios tend to similar 
values as for the hard-sphere and the LJ system. The 
ratios /tioo/mho and ^ioo/l*iii for the WCA model de¬ 
crease respectively from high values of 1.77 ± 0.099 and 
2.68 ± 0.13 at low pressures and then saturate at val¬ 
ues slightly smaller than those for the hard-sphere and 
the LJ system. These results indicate that the specific 
nature of the intermolecular potential has a substantial 
effect on the magnitude and anisotropy of /i and cannot 
be ignored, even though entropic effects play a dominant 
role. 

Now, we address the question to what extent clas¬ 
sical models of crystal growth can predict the depen¬ 
dence of the interfacial velocity on undercooling. The 



FIG. 10. (Color online) The kinetic growth coefficient fi cor¬ 
responding to the WCA potential for the three crystal orien¬ 
tations (100), (110), and (111), (a) as a function of the coex¬ 
istence pressure Pm, and (b) as a function of the coexistence 
temperature Tm- Note that in (b), jj,* = /t-v/Tm is shown. 
For clarity, in (a), /r for the (110) orientation is shown in the 
inset. Representation of symbols is same as in Fig. 


WF model |, ^ assumes that crystal growth is an ac¬ 
tivated process which is limited by the self-diffusion of 
the atoms in the liquid phase. This model leads to the 
following expression for EES: 

^;.WF ^ :^e(-AS/fcB) [1 _ e-AG/teT] (g) 

where, D is the self-diffusion constant of the liquid atoms, 
d the interplanar spacing between adjacent crystalline 
layers, A the mean free path of a liquid atom, / the prob¬ 
ability of a liquid atom to be attached to a crystal lattice 
site at the interface, and AS and AG respectively the 
differences in entropy and Gibbs free energy per particle 
between the crystal and liquid phases. The temperature 
is assumed to be below the melting temperature, T < Tm. 
The thermodynamic driving force for the crystallization 
is described by the term [1 — 

The WF model has been shown to predict the crystal¬ 
lization rates of binary (metallic) systems fairly well [ll|, 
Id^ . For one-component systems, however, the WF 
model tends to underestimate kinetic growth coeffi¬ 
cients. Therefore, Broughton, Gilmer and Jackson (BGJ) 
et al. 00 predicted an alternative collision-limited 
model, 


Vi 


BGJ 


# /3fc^e(-AS/fcB)ri _ g-AG/fcBTi 

A V m 


( 7 ) 


The BGJ model differs from the WF model in that the 
prefactor depending on the self-diffusion constant of the 
liquid atoms is replaced by a term containing their ther¬ 
mal velocity such that in the BGJ model the limiting fac¬ 
tor for crystal growth is the thermal velocity with which 
the atoms collide rather than the self-diffusion coefficient. 
Various works have compared these two theories to simu¬ 
lation and experimental results and it has been observed 
that the collision-limited model is an accurate predictor 
of fi (subject to an appropriate value for the fit param¬ 
eter /), at least for the (100) and (110) orientation of 
one-component fee systems 0, 0, 0 choosing the fit 
parameters 0 | / = 0.27 and A = 0.15a (with a the 
lattice constant). In simulations of LJ systems 00, 
however, the WF model has been shown to describe the 
interface velocity corresponding to the (111) orientation 
reasonably well. 

At small undercoolings, AG is proportional to the un¬ 
dercooling, AT, and the entropy difference can be taken 
as independent of temperature. Hence, close to coexis¬ 
tence 1 — exp(—AG)/A: bT Ri ksTTM ^ L/Tm- 

Inserting these expressions in Eqs. (jH) and ([7|) and us¬ 
ing the definition of /i, Eq. (ED, one obtains the following 
expressions for /i corresponding to the WF model, 

^WF ^ (g) 

and the BGJ model, 


A V m 


( 9 ) 
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FIG. 11. (Color online) The ratios /rioo//riio (circles) and 
Mioo/miii (squares) for the fsLJ system. The inset shows the 
same ratios for the WCA system. 
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FIG. 12. (Color online) The Lindemann ratio (in the inset, 
the RMS freezing ratio) as a function of pressure for the WCA 
(red squares) and fsLJ potentials (black circles). 


Note that since we consider small undercoolings, we can 
replace the temperature T in Eqs. dH]) and (0) by the 
melting temperature Tm. 

Now, we compare the above predictions for fj, with the 
simulation data along the P — T coexistence line. For 
both the fsLJ and the WCA model, A in Eq. is cho¬ 
sen to be 0.14a for comparison with results correspond¬ 
ing to the (100) and (110) crystal-liquid interfaces. The 
parameter / is chosen to be 0.35 (0.28) and 0.30 (0.29) 
for the (100) and (110) orientations of the fsLJ (WCA) 
model, respectively. To compare the WE model, Eq. 
to the simulation results for the (111) orientation of the 
fsLJ and the WCA model, the parameters A and / are 
chosen to be 0.04a and 0.27, respectively. Note that the 
same values of these parameters are used for each point 
along the coexistence line. 

Figures [5^ and [TUh show reasonable agreement be¬ 
tween theory and simulation. Using the classical mod¬ 
els, the decrease of the kinetic growth coefficient with 
increasing melting temperature can be understood in 
terms of the decrease of the thermodynamic driving 
force along the coexistence line, described by the term 
e(-^'S'/fcj3)j]^_g-AG/fcBTj Rj fcflTTM • The decrease 

of this term with increasing Tm overpowers the increase 
of the remaining terms containing D and i/Tm leading 
to slower growth kinetics. 

It is to be noted that the values of / in Eq. , cor¬ 
responding to the (100) and (110) orientations of the 
two potentials that leads to the best agreement with the 
simulation data, are very close, indicating that this pa¬ 
rameter is independent of the orientation. In Eq. (® [I 
corresponding to the different orientations are propor¬ 
tional to the interplanar spacing d yielding the value 
fiioo/Miio = 1.414. As shown in Fig. [TTl results from 
simulations are in reasonable agreement with this value. 


F. Freezing and melting rules 


The FSM method yields simultaneously the P — T 
coexistence line and the kinetic growth coefficient, /r. 
However, computational app roaches such as the capillary 
fluctuation metho d Il5l420l | or the equilibrium fluctua¬ 
tion technique [22l - l24| require an accurate determination 
of the coexistence conditions for obtaining precise esti¬ 
mates of the kinetic growth coefficient. Computing the 
coexistence conditions at several state points to obtain an 
accurate phase diagram, is computationally demanding. 
To reduce the computational effort, it would be helpful 
to have a prior idea of the crystal-liquid phase coexis¬ 
tence region, to avoid carrying out fruitless simulations 
at regions far away from the phase boundaries. To this 
end, there have been a number of attempts to identify 
certain empirical rules obeyed by the individual phases 
at melting and freezing [^. How accurately such 
rules are obeyed depends on the nature of the potential 
and the coexistence conditions. Here, we test the valid¬ 
ity of three empirical rules for the two models studied 
in this work: (i) Lindemann’s melting rule [s^, (ii) the 
Raveche-Mountain-Streett freezing rule [3^ and (iii) the 
Hansen-Verlet freezing criterion |40l |. 

According to Lindemann’s melting rule, a crystal melts 
when the root mean square displacement of crystalline 
atoms around their ideal lattice positions is approxi¬ 
mately 10% of the nearest-neighbour distance, a. To test 
this rule, we compute the Lindemann ratio L^, given by 


L, 


a 


( 10 ) 


The values of Tr, as reported from previous studies for 
different model potentials, range roughly from 0.1 to 0.2 

fU 


To compute Tr and also the quantities involved in the 
other empirical rules, we have carried out simulations of 





























10 


3.1 


3 


t 2.9 


2.8 


2.7 



WCAl 
* •fsU I 



10 


-2 


16^ 



10 


1(f 1(f 


FIG. 13. (Color online) The maximum of the liquid struc¬ 
ture factor, S'max(fc), along the coexistence line for the fsLJ 
(circles) and WCA (squares) system. 


The Hansen-Verlet freezing rule states that the first 
peak of the liquid structure factor attains the value 
<S'max(^) = 2.85 and remains invariant along the freez¬ 
ing curve. Hansen and Verlet postulated this rule based 
on observations corresponding to the LJ potential. Later, 
Agrawal and Kofke reported a 10% increase of this 
quantity along the coexistence line for the LJ potential. 
We computed the structure factor directly from the coor¬ 
dinates of the particles and not via a Fourier transform of 
the pair correlation function Q. Our results (Fig.[T31) for 
the two systems indicate a systematic increase of S'max(fc) 
along the freezing curve. While S'max corresponding to 
the fsLJ system increases from a minimum of 2.74 to 3.02 
at the highest pressure, that for the WCA potential re¬ 
mains above the value predicted by the Hansen-Verlet 
freezing criterion at all the considered coexistence pres¬ 
sures. This indicates that the Hansen-Verlet criteria can 
at most be used as a rule of thumb principle to indicate 
how close the system is to the actual freezing tempera¬ 
ture. 


the bulk crystal and liquid at the appropriate coexistence 
temperature and pressure for systems consisting of 8788 
and 6912 particles for the WCA and the fsLJ model, 
respectively. 

Figure [T^ shows the variation of along the coexis¬ 
tence line. For both models, it increases from a value 
of about 0.10 at low pressure to one of about 0.12 at 
high pressure, indicating that is not invariant along 
the melting line, but the closeness of the reported val¬ 
ues with Lindemann’s prediction support a qualitative 
validity of this rule. 

The second empirical rule we test is based on the liquid 
structure and known as the Raveche-Mountain-Streett 
(RMS) freezing criterion. According to this rule, along 
the freezing line, the radial distribution function of the 
liquid obeys the following relation: 

I = ff(?'inin)/3(?’max) = 0.2 ± 0.02 (11) 

where, girmin) is the first non-zero minimum and g(rmax) 
the first maximum of the radial distribution function 
g{r). In the inset of Fig. [1^ we show the RMS freez¬ 
ing criterion along the coexistence pressure for the WCA 
and the fsLJ model. While quantitatively our results dif¬ 
fer from the RMS rule, the data corresponding to the 
fsLJ model, which hover around 0.185 at low pressures 
(between 0.005 and I.O) and around 0.180 at larger val¬ 
ues of Pm, point to a rough invariance of this quantity 
along the coexistence line. However, at low pressures 
(Pm < 30), I corresponding to the WCA potential is 
much smaller than suggested by the RMS criterion, while 
at larger pressures, it saturates around value of about 
0.186, close to the value obtained for the fsLJ model. A 
similar variation of the RMS freezing criterion along the 
coexistence line was also observed in an earlier work . 


VI. CONCLUSION 

We have computed the crystal-melt kinetic growth co¬ 
efficient from molecular dynamics simulation via the free 
solidification method. Our results are consistent with 
previous works as regards the magnitude and anisotropy 
of g is concerned. The variation of the kinetic growth 
coefficient along coexistence indicates that the slowing 
down of the crystal growth is primarily a temperature 
effect. Changing the pressure by two orders of magni¬ 
tude in case of the fsLJ potential does not change the 
magnitude of the kinetic growth coefficient in the coex¬ 
istence region where the melting temperature remains 
almost unchanged. 

Similarity of values between Lennard-Jones and hard- 
sphere systems indicate that packing effects play a domi¬ 
nant role in describing the growth kinetics. However, the 
specific nature of the interaction potential cannot be ig¬ 
nored. Classical theories of crystallization models predict 
values of g in reasonably good agreement with simulation 
results, subject to an appropriate value of a fit parame¬ 
ter. In the high-temperature and high-pressure limit, the 
crystal growth kinetics of the fsLJ and WCA potentials 
are similar to each other. In this limit, the coexistence 
properties of the two systems approach that of the purely 
repulsive, inverse twelfth-power potential. Various melt¬ 
ing and freezing rules have been investigated and while 
they roughly indicate the coexistence region, they are not 
quantitatively accurate. 
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